/*
 * CCEengine.cpp
 *
 *  Created on: May 25, 2011
 *      Author: nzhao
 */
#include "share/Usrdef.h"

#include "core/globalvars.h"
#include "share/parameters.h"
#include "share/pulseseq.h"
#include "special/import_data.h"

#ifdef NV_CENTER
#include "special/material_NV.h"
#include "special/interaction_NV.h"
#endif

#ifdef SiC_4H
#include "special/material_4HSiC.h"
#include "special/interaction_4HSiC.h"
#endif

#include "share/GenClst_Control.h"
#include "share/HamiForm.h"
#include "core/Spin.h"
#include "core/Cluster.h"

char	WORKSHOP[200];
char	WORKSHOP_LOCAL[]	= "/home/nzhao/research/computation/workspace/cce-engine/";
char	WORKSHOP_SERVER[]	= "/home/user/nzhao/Projects/cce-engine/";
char	DEBUG_PATH[]		= "dbg/";

int		NBATH;						//number of bath spins
double	* COORDINATES;          	//bath spin coordinates
int		* ATOM_SPIN;
double	* ABSOLUTE_GAMMA_N;
int		BUFFER_CLST_NUM = 5000000;	//smaller than 5.0*10^6 clusters
double	* RES_BUFFER_REDUCE;
double	* RES_BUFFER_IRREDUCE;
//////////////////////////////////////////////////////////////////////////////////////////////////////////
/**\defgroup MAIN Main program
 * @{ */

struct ClusterInfo
{
	int nspin;
	double max_deviation;
	double sum_deviation;
};

clock_t		START,				/**< Program starting time*/
			FINISH,				/**< Program finishing time*/
			ABS_TIME_START,		/**< Absolute starting time*/
			ABS_TIME_END;		/**< Absolute finishing time*/
double 		TIME_USED,			/**< CPU time used*/
			ABS_TIME_USED;		/**< Absolute time used*/

double		* FINAL_SIGNAL;
double		** ORDER_RES;
int			PROGRESS;
ClusterInfo	* clstInfo;
cSPIN 		* BathSpinList;
cCLUSTER	CLUSTER;


void    Run_Time(const char * flag); /**< An assistent function for recording the time. */
void	CCE_Preparation();
void 	post_calculations(int i, int j, int count, double * clst_res, double * irreducible);
int 	export_RES(int argc, char **argv);
void 	collect_CCE(int argc, char **argv);
/**@}*/

int main(int argc, char **argv)
{
	int i, j, k, count;
    Run_Time("start");

    //set WORKSHOP path
	if (strcmp(WorkshopPath, "LOCAL") == 0)
		strcpy(WORKSHOP, WORKSHOP_LOCAL);
	else
		strcpy(WORKSHOP, WORKSHOP_SERVER);
	cout << "Workshop path = " << WORKSHOP << endl;

	ParameterResolve(argc, argv);

	set_time_sequence();

	import_data();

	calc_Interactions();

	generate_clusters();

    if(TASK_MODE == 1)
    	generate_initial_state(STATE_SEED);

    generate_hamiltonian_forms();

    /////////////////////////////////////////////////////////////////////////////////
    //CCE preparations

    BathSpinList = new cSPIN [NBATH];
    SingleSpinInitialization(BathSpinList);

    CCE_Preparation();

    /////////////////////////////////////////////////////////////////////////////////
    //cluster-by-cluster calculations

    double 	* clst_res = new double [NTIME + 1];
    double 	* irr_clst_res = new double [NTIME + 1];

    count = 0;
    for(i = 0; i < CALCULATED_ORDER; i ++)
    {
    	cout << "Calculating order " << (i + 1) << endl;
    	for(j = 0; j < NUM_CLST[i]; j ++)
    	{
    		if (count % PROGRESS == 0) cout << count / PROGRESS << "% completed." << endl ;

    		calculation_cluster(CLUSTER, BathSpinList, i, j, count, clst_res, irr_clst_res);
    		post_calculations(i, j, count, clst_res, irr_clst_res);

    		count ++;
    	}
		for (k = 0; k <= NTIME; k++) {ORDER_RES[i][k] = FINAL_SIGNAL[k];}
		cout << " completed! " << endl;
    }
    assert(count == TOTALCLST);

    /////////////////////////////////////////////////////////////////////////////////
    //finalize
    Run_Time("end");
    export_RES(argc, argv);
//    collect_CCE(argc, argv);
	return 0;
}
void CCE_Preparation()
{
    int     i;
    //////////////////////////////////////////////
    //allocate space: buffer
    int buffer_size, cache_size;
    if ( CACHE_Q )
    {
        cache_size = TOTALCLST;
        buffer_size = BUFFER_CLST_NUM > cache_size ? cache_size * (NTIME+1) : BUFFER_CLST_NUM * (NTIME+1);
        RES_BUFFER_REDUCE   = new double [buffer_size];
        RES_BUFFER_IRREDUCE = new double [buffer_size];
    }
    else
    {
        cache_size = TOTALCLST - NUM_CLST[CALCULATED_ORDER - 1];
        buffer_size = BUFFER_CLST_NUM > cache_size ? cache_size * (NTIME+1) : BUFFER_CLST_NUM * (NTIME+1);
        RES_BUFFER_IRREDUCE = new double [buffer_size];
    }
    cout << "CACHE_Q = "<< CACHE_Q << "; Cached cluster number = " << cache_size << endl;

    clstInfo = new ClusterInfo [TOTALCLST];
    for(i = 0; i < TOTALCLST; i ++)
    {
    	clstInfo[i].max_deviation = 0.0;
    	clstInfo[i].sum_deviation = 0.0;
    }

    //////////////////////////////////////////////
    //allocate space: results
    FINAL_SIGNAL= new double [NTIME+1];
    ORDER_RES 	= new double * [CALCULATED_ORDER];
    for (i = 0; i < CALCULATED_ORDER; i ++)
    {
    	ORDER_RES[i] = new double [NTIME+1];
    }

    switch(TASK_MODE)
    {
    	case 0:
    		/////////////////////////////////////////////////////////////////////////////////////
    		//Ensemble CCE
    		for(i = 0; i <= NTIME; i ++)
    			FINAL_SIGNAL[i] = 1.0;
    		break;
    	case 1:
			/////////////////////////////////////////////////////////////////////////////////////
			//Single-Sample CCE
    		for(i = 0; i <= NTIME; i ++)
    			FINAL_SIGNAL[i] = 1.0;
    		break;
    	case 2:
			/////////////////////////////////////////////////////////////////////////////////////
			//Correlation function
    		for(i = 0; i <= NTIME; i ++)
    			FINAL_SIGNAL[i] = 0.0;
    		break;
    	case 3:
			/////////////////////////////////////////////////////////////////////////////////////
			//Gaussian approximation
    		for(i = 0; i <= NTIME; i ++)
    			FINAL_SIGNAL[i] = 1.0;
    		break;
    	case 4:
			/////////////////////////////////////////////////////////////////////////////////////
			//Noise spectrum
    		for(i = 0; i <= NTIME; i ++)
    			FINAL_SIGNAL[i] = 0.0;
    		break;
    	default:
			/////////////////////////////////////////////////////////////////////////////////////
			// error
			cout << "ERROR: Wrong TASK_MODE" << endl;
			assert(0);
			break;
    }

    //////////////////////////////////////////////
    //progress

    PROGRESS = TOTALCLST / 100 + 1;

}

void post_calculations(int order, int idx_j, int clst_idx, double * reducible, double * irreducible)
{
	double diff, max_diff, sum_diff;
	int i;

	switch (TASK_MODE) {
		case 0:
			/////////////////////////////////////////////////////////////////////////////////////
			//Ensemble CCE
			max_diff = 0.0;	sum_diff = 0.0;
			for(i = 0; i <= NTIME; i++)
			{
				FINAL_SIGNAL[i] *= irreducible[i];

				diff = fabs( irreducible[i] - 1.0);
				sum_diff += diff;
				max_diff = (diff > max_diff ? diff : max_diff);
			}
			clstInfo[clst_idx].max_deviation = max_diff;
			clstInfo[clst_idx].sum_deviation = sum_diff;

			break;

		case 1:
			/////////////////////////////////////////////////////////////////////////////////////
			//Single-Sample CCE
			max_diff = 0.0;	sum_diff = 0.0;
			for(i = 0; i <= NTIME; i++)
			{
				FINAL_SIGNAL[i] *= irreducible[i];

				diff = fabs( irreducible[i] - 1.0);
				sum_diff += diff;
				max_diff = (diff > max_diff ? diff : max_diff);
			}
			clstInfo[clst_idx].max_deviation = max_diff;
			clstInfo[clst_idx].sum_deviation = sum_diff;

			break;

		case 2:
			/////////////////////////////////////////////////////////////////////////////////////
			//Correlation function
			max_diff = 0.0;	sum_diff = 0.0;
			for(i = 0; i <= NTIME; i++)
			{
				FINAL_SIGNAL[i] += irreducible[i];

				diff = fabs( irreducible[i] - 0.0);
				sum_diff += diff;
				max_diff = (diff > max_diff ? diff : max_diff);
			}
			clstInfo[clst_idx].max_deviation = max_diff;
			clstInfo[clst_idx].sum_deviation = sum_diff;

			break;

		case 3:
			/////////////////////////////////////////////////////////////////////////////////////
			//Gaussian approximation
			max_diff = 0.0;	sum_diff = 0.0;
			for(i = 0; i <= NTIME; i++)
			{
				FINAL_SIGNAL[i] *= irreducible[i];

				diff = fabs( irreducible[i] - 1.0);
				sum_diff += diff;
				max_diff = (diff > max_diff ? diff : max_diff);
			}
			clstInfo[clst_idx].max_deviation = max_diff;
			clstInfo[clst_idx].sum_deviation = sum_diff;

			break;

		case 4:
			/////////////////////////////////////////////////////////////////////////////////////
			//Noise spectrum
			max_diff = 0.0;	sum_diff = 0.0;
			for(i = 0; i <= NTIME; i++)
			{
				FINAL_SIGNAL[i] += irreducible[i];

				diff = fabs( irreducible[i] - 0.0);
				sum_diff += diff;
				max_diff = (diff > max_diff ? diff : max_diff);
			}
			clstInfo[clst_idx].max_deviation = max_diff;
			clstInfo[clst_idx].sum_deviation = sum_diff;

			break;

		default:
			/////////////////////////////////////////////////////////////////////////////////////
			// error
			cout << "ERROR: Wrong TASK_MODE" << endl;
			assert(0);
			break;
	}

}

int export_RES(int argc, char **argv)
{
    int i, j;
    double tval;

    ///////////////////////////////////////////////////////////////////////////////////
    //print to screen
    for (i = 0; i <= NTIME; i++)
    {
        cout << MIN_TIME + i * TIME_STEP  << "\t"
            << setw(12) << setiosflags(ios::scientific) << setprecision(10)
            << FINAL_SIGNAL[i] << endl;
    }

    ///////////////////////////////////////////////////////////////////////////////////
    //save to file
    char res_file[200];

    strcpy(res_file, OUTPUT_PATH);
    strcat(res_file, OUTPUT_FILE_NAME);
    strcat(res_file, ".dat");

    ofstream res(res_file);
    if(!res) assert(0);

    res << "#info\t" << "Start Time: " << ctime(&ABS_TIME_START) ;
    res << "#info\t" << "Finish Time: " << ctime(&ABS_TIME_END) ;
    res << "#info\t" << "Absolute Time: " << ABS_TIME_USED << " seconds." << endl;
    res << "#info\t" << "CPU Time: " << TIME_USED << " seconds." << endl;
    res << "#info\t" << "VERSION: " << VERSION_INFO << endl ;
    res << "#info\t" << "GROUP ID: " << GROUP_ID << endl;
    res << "#info\t" << "Exec: " << EXEC << endl;

    res << "#para\t" << "Input file:\t"	    << INPUT_FILE_NAME << endl;
    res << "#para\t" << "Config. file:\t"   << CONFIG_FILE_NAME << endl;


    res << "#para\t" << "Extern magB:\t"     << EXT_MAG_FIELD << endl;
    res << "#para\t" << "Theta angle:\t"     << THETA << endl;
    res << "#para\t" << "Phi angle:\t"       << PHI << endl;
    res << "#para\t" << "Pulse type:\t"       << STR_SEQ_NAME << endl;
    res << "#para\t" << "Pulse number:\t"     << PULSE_NUM << endl;
    res << "#para\t" << "Is RF Driving:\t"	<< IsDriving << endl;
    if(IsDriving)
    {
    	res << "#para\t" << "Drv Detuning:\t"	<< DrvDetuning << endl;
    	res << "#para\t" << "Drv Rabi:\t"		<< DrvRabi<< endl;
    }
    res << "#para\t" << "CCE order:\t"        << CALCULATED_ORDER << endl;
    for (i = 0; i < CALCULATED_ORDER; i++)
    {
        res << "#para\t" << "Cluster Number of Order " << i + 1 << ":\t" << NUM_CLST[i] << endl;
    }
    res << "#para\t" << "Time start:\t"      << MIN_TIME << endl;
    res << "#para\t" << "Time end:\t"        << MAX_TIME << endl;
    res << "#para\t" << "Time point:\t"      << NTIME << endl;
    res << "#para\t" << "Threshold:\t"       << sqrt(THRESHOLDVALUE)<< endl;
    res << "#para\t" << "Transition:\t"      << STATE_1 << "--" << STATE_2 << endl;
    res << "#para\t" << "CACHE_Q:\t"	     << CACHE_Q << endl;
    res << "#para\t" << "Bath Spin Num:\t"   << NBATH << endl;
    res << "#para\t" << "IsCore:\t"			 << IsCore << endl;
    res << "#para\t" << "IsBlock:\t"		 << IS_BLK << endl;

    res << "#axes\t" << "x\t";
    for (i = 0; i < CALCULATED_ORDER; i ++) res << "y\t";
    res << endl;

    res << "#name\t" << "time\t";
    for (i = 0; i < CALCULATED_ORDER; i ++) res << "signal\t";
    res << endl;

    res << "#unit\t"<< "ms" << "\t" ;
    for (i = 0; i < CALCULATED_ORDER; i ++) res << "a.u.\t";
    res << endl;

    res << "#lgnd\t" << "t\t" ;
    for (i = 0; i < CALCULATED_ORDER; i ++)
    {
        res << STR_SEQ_NAME << PULSE_NUM <<"order" << i+ 1;
        if(i == CALCULATED_ORDER - 1)
        {res << endl;}
        else
        { res << "\t";}
    }

    for (i = 0; i <= NTIME; i++)
    {
        tval = MIN_TIME + i * TIME_STEP;
        if(TASK_MODE == 4)
        	tval /= (2.0 * PI);
        res << setw(12) << setiosflags(ios::scientific) << setprecision(10);
        res << tval << "\t";
        for (j = 0; j < CALCULATED_ORDER; j++)
        {

            res << ORDER_RES[j][i] ;
            if(j == CALCULATED_ORDER - 1)
            {res << endl;}
            else
            { res << "\t";}
        }

    }

    res.close();

    if ( CACHE_Q )
    {
    	char cluster_data_file[200];
    	strcpy(cluster_data_file, OUTPUT_PATH);
    	strcat(cluster_data_file, "Cluster_");
    	strcat(cluster_data_file, OUTPUT_FILE_NAME);
    	strcat(cluster_data_file, ".dat");

    	ofstream clstdata(cluster_data_file, ios::binary);
    	if(!clstdata) assert(0);

    	clstdata.write((const char *)&TOTALCLST, sizeof(int));
    	clstdata.write((const char *)&CALCULATED_ORDER, sizeof(int));
    	clstdata.write((const char *)&NTIME , sizeof(int));
    	clstdata.write((const char *)NUM_CLST , CALCULATED_ORDER * sizeof(int));
    	clstdata.write((const char *)RES_BUFFER_IRREDUCE, TOTALCLST * (NTIME+1) * sizeof(double) );

    	clstdata.close();
    }


    return 0;
}

void collect_CCE(int argc, char **argv)
{
	int i;
	cout << "Calculating the cluster amplitude......" ;

	char singular_file[200];
	strcpy(singular_file, OUTPUT_PATH);
	strcat(singular_file, "singular");
	strcat(singular_file, ".dat");

	ofstream singular(singular_file);
    if(!singular) assert(0);

    for(i = 0; i < TOTALCLST; i++)
    {
    	singular << i << "\t" << clstInfo[i].max_deviation << "\t" << clstInfo[i].sum_deviation << endl;
    }


	singular.close();

}

void Run_Time(const char * flag)//(int k)
{
    //if(k==1)    //start
    if(strcmp(flag, "start") == 0 )
    {
        START = clock();
        time(&ABS_TIME_START);
    }
    else if(strcmp(flag, "end") == 0 )
    {
       	FINISH = clock();
       	time(&ABS_TIME_END);
        TIME_USED = (double)(FINISH - START) / CLOCKS_PER_SEC;
        ABS_TIME_USED = difftime(ABS_TIME_END, ABS_TIME_START);
        cout << "calculation time:" << endl;
        cout << "CPU time = " <<  TIME_USED << " second;" << endl;
        cout << "Absolute Time = " << ABS_TIME_USED << " second." <<   endl;
    }
    else
    {
        cout << "do nothing!" << endl;
    }
}

///////////////////////////////////////////////////////////////////////////////////////
/*!\fn void Run_Time(const char * flag)
 * \param flag A string, either "start" or "end", which record the start or end time.
 * Other input parameters will be ignored.
 *
 * This function records the program starting or ending time.
 *
 * If \a flag=="start", \a START, and \a ABS_TIME_START will be set.
 *
 * If \a flag=="end", \a FINISH, \a ABS_TIME_END will be set, and cup time and absolute time is calculated and printed.
 */

///////////////////////////////////////////////////////////////////////////////////////
/*!\var clock_t E
 * The program starting time, obtained by \a clock().
 */

///////////////////////////////////////////////////////////////////////////////////////
/*!\var clock_t FINISH
 * The program finishing time, obtained by \a clock().
 */

///////////////////////////////////////////////////////////////////////////////////////
/*!\var clock_t ABS_TIME_START
 * The Absolute program starting time, obtained by \a time(&ABS_TIME_START).
 */

///////////////////////////////////////////////////////////////////////////////////////
/*!\var clock_t ABS_TIME_END
 * The Absolute program finishing time, obtained by \a time(&ABS_TIME_END).
 */

///////////////////////////////////////////////////////////////////////////////////////
/*!\var double TIME_USED
 * The CPU time of the program in second.
 */

///////////////////////////////////////////////////////////////////////////////////////
/*!\var double ABS_TIME_USED
 * The absolute time of the program in second.
 */
